Epidemic threshold and control in a dynamic network 
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In this paper we present a model describing Susceptible-Infected-Susceptible (SIS) type epidemics 
spreading on a dynamic contact network with random hnk activation and deletion where link ac- 
tivation can be locally constrained. We use and adapt a improved effective degree compartmental 
modelling hamework recently proposed by Lindquist et al. [J. Lindquist et ai, J. Math Biol. 62, 2, 
143 (2010)] and Marceau et al. [V. Marceau et al, Phys. Rev. E 82, 036116 (2010)]. The resulting 
set of ordinary differential equations (ODEs) is solved numerically and results are compared to those 
obtained using individual-based stochastic network simulation. We show that the ODEs display ex- 
cellent agreement with simulation for the evolution of both the disease and the network, and is able 
to accurately capture the epidemic threshold for a wide range of parameters. We also present an 
analytical _Ro calculation for the dynamic network model and show that depending on the relative 
timescales of the network evolution and disease transmission two limiting cases are recovered: (i) 
the static network case when network evolution is slow and (ii) homogeneous random mixing when 
the network evolution is rapid. We also use our threshold calculation to highlight the dangers of 
relying on local stability analysis when predicting epidemic outbreaks on evolving networks. 
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I. INTRODUCTION 

The rise in the popularity and relevance of networks 
as a tool for modelling complex systems is well illus- 
trated by the ever increasing body of research concerned 
with the spread of diseases within host populations ex- 
hibiting non-trivial contact structures [l|, |2|. Networks 
offer an intuitive and relatively simple modelling frame- 
work which enables us to relax the strong implicit as- 
sumptions of more classical ordinary differential equa- 
tions (ODE) based approaches and to account for com- 
plexities in the contact structure of the host population 
[3-0] ■ This approach has shown that epidemic thresholds 
not only depend upon the infectiousness of the pathogen, 
or even simply the mean number of contacts per individ- 
ual, but also upon the exact structure of the host popu- 
lation [1,1^. In addition to its inherent theoretical value, 
this paradigm has immediate practical benefits, as the 
primary role of public health services is to put measures 
in place to bring diseases below their epidemic thresh- 
old. These measures depend heavily upon disrupting the 
transmission of a disease through vaccination and also 
more directly through the closure of public services, or 
even quarantine and curfews in extreme cases. Hence 
the knowledge of how the structure of the host popula- 
tion is contributing to the spread of a disease would help 
to increase the efficacy of any intervention [To| . 

Despite advances in both rigorous and non-rigorous 
analysis of networks, a key assumption in many network 
models is that contacts are fixed for the duration of an 
epidemic and that the disease propagates with a constant 
intensity across links. This will not be true for many dis- 
eases, especially those with long infectious periods, or 
diseases that become endemic. Indeed human contact 
patterns are well described by short repeated events, with 
individuals having a number of contacts best described by 



some appropriate time dependent random variable 
Furthermore, individuals and the communities they be- 
long to are likely to change their contact behaviour as a 
result of natural evolution and endogenous or exogenous 
perturbations such as a disease outbreak [l^ . 

Recently a number of studies have attempted to re- 
lax this assumption by allowing the networks to evolve 
over time by either v ary i ng contacts independently of the 
status of individuals [13, 13 or by explicitly coupling con- 
tact activation and deletion to the disease status of in- 
dividuals Thus, in the latter case, the dynamics 
of the disease is coupled with the dynamics of the net- 
work itself, with both acting as a feedback mechanism for 
the other [16,, il^, [l^ . Many of these studies have built 
macro ODE-based models that describe the coevolution 
of networks and the diseases that spread along them [Tsl- 
fl?. . 20] . All these studies confirm that dynamic networks 
and the coupling between the two dynamics lead to a 
richer spectrum of behavior than is found for epidemics 
on static networks. 

A crucial feature of allowing the co-evolution of dis- 
ease and network is the interplay and feedback between 
both dynamics, however this interdependence is diffi- 
cult to measure empirically. The models developed so 
far mainly use rewiring rules that intuitively make sense 
given that individuals would have knowledge of the dis- 
ease states of the rest of the population. However in this 
paper we move away from these assumptions and we pro- 
pose a dynamic network model that is based on random 
link activation-deletion, which would be more relevant 
for asymptomatic diseases, such as Chlamydia i21|]. Fur- 
thermore our dynamic network model is refined by intro- 
ducing a local constraint on link activation to account for 
the difference in the magnitude of the number of contacts 
of a node relative to system size. This dynamic network 
coupled with the simple Susceptible-Infected-Susceptible 
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(SIS) disease dynamics leads to the full model that will 
be analysed and discussed. We study this system and ex- 
plore to what extent a macro ODE-based compartmental 
model proposed for static networks is flexible enough to 
be adapted to a dynamic network case. Specifically, we 
focus on the SIS effective degree model as described in 
detail by Lindquist et al. [22| and also, to our knowl- 
edge, proposed by Marceau et al. [l^ in close succession. 
Gleeson [2j| later uses this same modelling framework 
and demonstrates that the effective degree formulation 
can be used to model other binary-state dynamics such 
as Glauber spin dynamics and shows that the ODE model 
can be used to carry out linear stability type analysis. 

Whereas both Lindquist et al. [13] and Gleeson [23| 
confine themselves to modelling on static contact net- 
works, Marceau et al. jT^] uses this same improved ef- 
fective degree formalism to explore SIS disease dynamics 
on adaptive networks. In this model the number of links 
in the network is fixed but the susceptible individuals 
can replace links to infectious neighbours with links to 
other randomly chosen susceptible individuals, as origi- 
nally proposed by Gross et al. . Our proposed model 
also uses SIS type epidemics on dynamic networks, but 
unlike Marceau our model allows for the random activa- 
tion and deletion of links over time. As such not only 
the network topology will evolve evolve and change over 
time, but also the number of links. This modified dy- 
namic effective degree model is also governed by a closed 
set of ODEs, which is then solved and compared to re- 
sults from individual based simulations and its ability to 
accurately predict the epidemic threshold over a range of 
parameters is investigated. We also derive an analytical 
Rq calculation that describes the stability of the disease- 
free equilibrium and we discuss the limitations of such 
a calculation in the light of having a dynamically active 
and evolving contact network. 



II. THE MODEL 

Linquist et al. [22| and Marceau et al. fl7j use different 
notation to describe the same modelling framework. For 
consistency, in this paper we follow the notation used by 
the former throughout. The effective degree modelling 
approach for SIS type disease dynamics not only 
categorizes the disease state of each individual as sus- 
ceptible {S) or infected (/) but also describes the state 
of their immediate neighbourhood. This is achieved by 
keeping track of the number of susceptible and infected 
neighbours that belongs to a given node. For example, 
Ssi represents the number of susceptible individuals that 
have s susceptible and i infected neighbours. This gives 
rise to more states and equations than would be seen in 
a standard pairwise model, where equations are given at 
the population level for all types of singles and pairs [13] . 
For example if a Sgi type node became infected via one 
of its i infectious neighbours, this individual would move 
to state Isi as only the status of the node itself is chang- 
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FIG. 1: (Colour online) Flow chart showing transitions in 
the dynamic SIS effective degree model. The directed red 
(gray), green (light gray), blue (dark gray) and black lines 
represent changes in state of an individual via infection, re- 
covery, link creation and link deletion respectively. The thick 
lines represent changes to the individual, and thin lines repre- 
sent changes to that individual's immediate neighbourhood. 
In relation to nodes of type Xsi,X £ {5,7}, infection of 
neighbours occurs at rate sGx, recovery of neighbours at 
rate 74, creation of a susceptible (infectious) link at rate 
a{AI — (s + i))Ps(i) and deletion of a susceptible (infectious) 
link at rate i^s(i), where: 



p ^k = l ^j + l = k 



, ,^ k=l 2^j + l = k ' 

'^i — P-^Tn — ?^ 77 — 



and = £f=oS, + i = fc(M-U+0)X„ 

Ef=o T., + i = kiM-(i+l))(.S^i + Iji) ■ 



ing. However, if one of the i infected neighbours of an 
Ssi type node recovered then the node would enter the 
Ss+i,i-i class, whereas infection of one of the s neigh- 
bouring susceptible nodes moves the Ssi type node into 
the Ss-i,i+i class. 

Lindquist et al. [l^l defined 7 to be the per node recov- 
ery rate, /3 the per link infection rate and M the max- 
imum nodal degree of a network with N nodes. They 
then derived the following system of ^^f^j^ 2(fc + 1) — 
M [M + 3) equations: 



Ssi = -l3iSsi+'-iIsi+l[{i + ^)S. 
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(1) 



^ n^k=lJ2j+l=kj^Sjl - , 



J2k=l ^j+l=k jSjl 

Isi = liiSsi - jisi + 7[(* + ^)Is~i,i+i - ill 



(2) 



r.^k = l^J+l = k^^Sjl 



for {(s, i) : s,i > 0,1 < k ^ s + i < M}. This is the SIS 
effective degree model for a static contact network. 
In oder to adapt this model to describe SIS dynamics 
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on a dynamic contact network, we introduce two new 
parameters: w, the per link deletion rate and a, the per 
non-link, or more precisely the per potential link creation 
rate. These rates could also be made to be link-type 
dependent, i.e. ljsi would be the per SI link deletion 
rate. For the dynamic network cae, the system size will 
increase slightly from M{M+3) to J2tLo 2(^+1) = {M+ 
l){M + 2) equations to account for nodes of the type Xq q 
where X E {S, I}. In the static case, these nodes were 
dynamically unimportant as they could neither infect nor 
become infected by other nodes. However in the dynamic 
model, they could connect to other nodes in the system 
and so enter states Xi^o or Xq^i depending on the state 
of the node with which they have just formed a new link. 

The total number of links in the system at time t, A{t), 
and potential links, ^{t) can easily be calculated from the 
effective degree formulation as 
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with the mean nodal degree given by {k{t)) = At 
the equilibrium, = loA which gives us the mean nodal 
degree: 



{ky 



a + u 



-M. 



(3) 



Note that Eq. does not depend on the system size, 
N, but rather on the maximum nodal degree, M. This is 
important because in the static model, M is simply given 
by the node or nodes with the highest degree whilst in 
the dynamic case, however, M can be considered as a 
carrying capacity, whereby no node can have more than 
M links. This subtle but important difference means 
that in the dynamic case, M itself can be regarded as a 
parameter which controls the potential level of network 
saturation. 

When adding the terms that govern link creation and 
deletion to Eqs. ([T]) and ([2]) it is far simpler to construct 
the terms that govern deletion of existing links than those 
for the creation of new links. Links to nodes of type X^i 
where X 6 {S, 1} are cut at a rate proportional to their 
degree, so individuals will leave Xgi through link dele- 
tion at a rate uj{s + i) and will either enter the Xs~i^i 
or Xs^i~i classes depending on the state of the nodes to 
which they were previously connected. Similarly individ- 
uals can enter state X^i if they were in states Xs^i+i or 
Xs+i,i and a link to an infected or susceptible node was 
deleted respectively. 

When creating new links to nodes of type X^i, there 
are M — {s + i) stubs remaining, so nodes will transi- 
tion out of this state at a rate a(M — (s -t- i)) and will 
either enter the Xs+is or Xs^i+i classes depending on 



FIG. 2: (Colour online) Time evolution of I{t) = 

H'k^oY^j+i^k^oii^) ^^'^ {k){t) = ^ for three different val- 
ues of M. Results from the ODE are given by solid lines 
and those from simulation by points. In all cases A'^ = 1000, 
Jo = 100, a = 0.05, CO = 0.1, P = 0.5 and 7=1. The initial 
network is a regular random graph with k = 4. In each case, 
mean values from the stochastic simulations were found by av- 
eraging over 100 repetitions, with the individual realisations 
plotted in grey. 



the state of the node to which they have just connected. 
The rate at which nodes enter the X^i class from ei- 
ther Xs-i^i or Xs^i-i depends not only on the number 
of stubs still available in the node in question, but also 
on the probability that the newly created link attaches 
to a node of state S or I respectively. So nodes enter 
Xsi from Xs-i.i at the rate aPs{M — {s ~ 1 + i)), and 
nodes enter Xsi from Xg^i-i at rate aPi{M — {s + i — 1)), 
EfcLo E,+,=JM-(3+o)x,, 



the probability of picking an available stub belonging to 
nodes of type X where X G {S, I}. The full set of tran- 
sitions captured by this model is shown in Fig. [TJ 

The addition of these terms to Eqs. (HI and ([2]) trans- 
forms the SIS effective degree model for a static network 
into one that captures the spread of SIS type diseases 
on a dynamic contact network and is described by the 
following system of {M + 1)(M + 2) equations: 



Ssi = -f3iSsi+'jIsi+l[{i + i)Ss-i^i+i~iSs 



(4) 



J2k=o J2j+i=k j^ji 

-uj[{s + i)Ssi - (i + - (s + l)S's+i,,;] 

-a{M - {s + i))Ss^ + a{M - (s - 1 + ^))PsSs^l,^, 
+a{M -{s + i-l))PISs,^-l 

l3iSsi - "fist + + ~ ilsi] (5) 
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-w[(s + i)Isi - {i + l)Is,i+l - (s + l)/s+i,,;] 
-a{M - (s + i))Is^ + a{M - (s - 1 + i))PsIs-i,i 

+a{M {s + i^l))PJIs,^-l, 

for {(s, i) : s,i > 0,0 < k — s + i < M}. This system is 
the dynamic SIS effective degree model. 



III. CALCULATING THE DISEASE 
THRESHOLD 

For the static case, Lindquist et al. [13] used the next 
generation matrix approach [IH to calculate the disease 
threshold to be 
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In this approach, Eqs. (|H) and (0 are linearized at the 
disease-free equilibrium (DFE) and the Jacobian at the 
DFE is written a.s F — V. In this formulation, F accounts 
for transitions from disease-free states to disease states 
(in the static case only the transition from Sgfl to 
needs to be considered) and V accounts for transitions 
between different disease states. The spectral radius, p, 
the leading eigenvalue of FV~^, gives Rq and describes 
the stability of the DFE. If i?o < 1 the DFE is stable 
and no epidemic will occur, but if i?o > 1 the DFE is 
unstable and the infectious agent can spread through the 
population. 

We can calculate F in the dynamic case by noting that 
the same 5^,0 to 5*^-1,1 type transitions can still occur, 
but in addition nodes can enter the disease states by link- 
ing to an infected node, namely S'^.o to Sg^i transitions. 
If we introduce a subscript s to denote the static version 
of the next generation matrix, so the static version of F 
is called Fg and so on, we have 



T T 

V V 



(7) 



where u^^. and Vg^. are {2k + 1) x 1 vectors. The u^^, vec- 
tors have kSk.o a-s their first entry and zeros elsewhere 
and the Vg^. vectors have their first (fc — 1) entries equal 
to {k — 1), 2(fc — 2), ... , s{k — s), . . . , (fc — 1) and zeros 
elsewhere. This is almost identical to the F matrix con- 
structed by Lindquist et al., but is augmented by Ugg 
and to account for the new disease state, /o,o, and 
the summation starts at fc = rather than fc = 1, 

We now introduce a new subscript d to describe the 
new transitions that are only possible in the dynamic 



FIG. 3: (Colour online) Epidemic threshold plot in the (/3, u) 
parameter space for four distinct values of a. Results from 
the ODE are given by solid lines and those from simulation 
by solid points. In each case, N = 1000, lo = 10, M — 20 and 
7 = 1. The initial network is a regular random graph with 
= 4. 



model. Hence a new F matrix, F^, is created, which has 
exactly the same dimensions as Fg, and is given by 



"do 



F, 



k)Skfi 



(8) 



Here, Ud^ is again a (2fc -|- 1) x 1 vector with the first 
entry equal to (M — (fc — l))5'fc_i^o and all other entries 
equal to zero. In the case where fc = 0, Udg = (0). In 
addition, Vd^. is the same size as Ud^. and the first fc entries 
are equal to zero, with the remaining fc+1 entries equal to 
Af — fc. The final F matrix that captures all the possible 
transitions in the dynamic effective degree model is found 
by taking a linear sum of the two, namely F = Fg + Fd- 

As with the static case, the V matrix is constructed 
through careful book-keeping, which can be done through 
iterative routines. In the static case, as the nodes have 
fixed degree, Vg is a block diagonal matrix with Vg — 
I4i ® © ■•• ffi VsM- For the dynamic model, Vd will be 
a block tri-diagonal matrix, as state transitions can now 
also occur by nodes gaining or loosing a link. In addition, 
the extra disease state /o,o needs to be considered, and 
V will now also depend upon a and w as well as /3 and 
7. Once F — Fg + Fd and V = Vd are constructed, the 
leading eigenvalue or Rq is computed numerically. 
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FIG. 4: (Colour online) Threshold stability in the -Ro) 
space with 7 = 1, Af = 20 and (fc)* — 3 for (thin solid 
lines, in order from top to bottom) a = 10~* (green), a = 
10"^ (red), a = 10"^ (blue) and a = 10 (black). In (a) the 
initial network is a regular random graph with fc = 6 and in 
(b) the initial degree distribution is negative binomial with 
(k) = 6 and = 12. In each case, cj = a^^^j^^. The thick 
short-dashed red line is the theoretical value of Rq for a static 
network, and the thick red dash-dotted line is the mean field 
limit 7?o = -(fe)*. 



IV. RESULTS AND DISCUSSION 

As shown in Fig. [21 the ODEs given by Eqs. (|4]) and 
closely capture the time evolution of an epidemic as 
predicted by stochastic simulations. The only parameter 
that is varied in Fig. [5] is M, and it is interesting to note 
the effect it has on the evolution of the disease. As per 
Eq. the mean nodal degree at equilibrium is depen- 
dent on M, and hence, given the same initial network 
configuration and values of a and uj, the network either 
loses or gains links as the system evolves. Thus varying 
the carrying capacity alone leads to different outcomes 
depending on whether the network can reach a level of 
connectedness that allows an epidemic to spread and be- 
come established. Allowing M to become an active model 
parameter that is able to control the outcome of an epi- 
demic has potentially interesting real world implications. 
The number of contacts per person is a natural, count- 
able property unlike the other model parameters, such as 
w, which are more difficult to infer. Therefore local con- 
straints that limit the maximum number of contacts per 
person could be potentially used as a metric when pro- 
moting safe behaviour at a population level in the event 
of an outbreak or other public health crisis. 

In Fig.[3l for a given value of a, M and /?, the epidemic 
threshold has been calculated from the ODEs in terms of 
id and compared to that predicted by simulations. The 
agreement is excellent and this is strong evidence that the 



FIG. 5: (Colour online) Time evolution of I and (fc) with 
7 = 1, Af = 20, a = 0.1, u = a '^''/^Jy} and an initial regular 
random network with fc = 6. The two cases illustrated above 
correspond to: (k)* — 3 and j3 = 0.35, giving Rg « 1.29 
(red long-dashed line) and (fc)* — 9 and /3 = 0.125, giving 
Rq ~ 0.77 (blue short-dashed line). 

dynamic effective degree model accurately captures the 
evolution of an epidemic on a network with random link 
creation and deletion. When considering the (/3, w) pa- 
rameter space used for the threshold plot in Fig. [3l there 
are three distinct regions that are worth noting. Firstly, 
given an initial starting network, it is possible to calcu- 
late the threshold value of /? in the static network case. 
For the regular random graph with fc = 4 used here, that 
value is (3* « 0.36. For values of /3 < 0.36, the relative 
time scales of disease and network evolution are crucial 
in determining whether or not an epidemic will occur. In 
this situation, the network needs to quickly evolve to be- 
come more densely connected in order for there to be an 
outbreak. The second area of interest is when the disease 
is highly infectious and as a result requires a high value 
of UJ to drive the epidemic below threshold. Indeed, if the 
disease parameters f3 and 7 are fixed then the only way of 
affecting the outcome of an epidemic is through changing 
the network structure, i.e. reducing the number of links 
or the variance. Hence, for a fixed a and M, a value of /3 
can be chosen large enough so that the minimum value of 
UJ needed to reduce the connectivity of the network suffi- 
ciently to stop an outbreak (see FigIS]), gives (fc)* < 2 as 
can be calculated from Eq. [31 If a network has (fc)* < 2 
then it becomes fragmented, with many nodes becoming 
unconnected. In these situations, the value uj needed to 
prevent an epidemic virtually destroys the network. In 
terms of real world implications, a large value of uj could 
correspond to a situation of strict quarantine and curfew 
whereby links between individuals are kept to a mini- 
mum. In between these two cases lies a region within 
which an epidemic would take hold naturally, given the 
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initial network, but which can be prevented by a value of 
Lu that leaves the network well connected. 

In Fig. m we show analytical values of i?o for a range 
of values of /3 and a. It is worth noting that two limiting 
cases are recovered when the timescale of the network dy- 
namics is fast and slow relative to the timescale of the dis- 
ease dynamics. The thick short-dashed red line shows Rq 
calculated for a static network, as proposed by Lindquist 
et al. [in and given in Eq. |6l and this is exactly followed 
by results from our dynamic Rq calculation when the net- 
work dynamics are set to be much slower than the disease 
dynamics. The other extreme is shown by the thick dash- 
dotted red line, and is the value of i?o that results form 
the classic mean- field calculation i?o — -^^p-- The time 

evolution of {k) is given by (fc) — a{M — (k)) —ijj{k) but, 
when the network dynamics is fast, the equilibrium net- 
work distribution, and hence (fc)*, is approached much 
quicker than the epidemic timescale and hence a value of 
(fc) = (fc)* as given by Eq. [3] can be used. This limit is 
closely matched by results from our dynamic Rq calcula- 
tion when the network dynamics are rapid compared to 
disease transmission as shown in Fig. 21 

Although Fig. |4] demonstrates the accuracy of our ana- 
lytical i?o calculation, Fig.[5]highlights two example cases 
where the long term epidemic outcomes are the oppo- 
site of what is predicted by i?o- In the cases i?o < 1 
(blue short-dashed curve) and i?o > 1 (red long-dashed 
curve) the system settles to an endemic and to a disease 
free equilibrium respectively, due to the different ways 
the networks evolve. Given that i?o is based on a local 
stability analysis, it can only incorporate the immedi- 
ate next-generation effects of random link activation and 
deletion, and cannot account for long term changes to 
the network structure. It is well established in the litera- 
ture (see, for example Li et al. '26']) that i?o is of limited 
value when used as a predictor, and even for static net- 



works needs to be used with care. Our results add weight 
to this argument, and we show that when dealing with 
disease spreading through dynamic contact networks the 
use of Ro as any kind of predictor on long term disease 
evolution should be met with some degree of caution. 

In summary, this paper has proposed an effective de- 
gree model for epidemics on dynamic networks with ran- 
dom link activation and deletion, where activation is lo- 
cally constrained. We have shown that this model agrees 
extremely well with results obtained from stochastic sim- 
ulations, and as such can reliably be used for the analyti- 
cal and semi-analytical study of coupled disease and net- 
work dynamics. We have shown how a local constraint 
limiting the number of contacts per individual can be 
used to control and prevent the outbreak of an epidemic 
in this dynamic model. We have also proposed an ana- 
lytical calculation of i?o, but also demonstrated the lim- 
ited value of threshold stability analysis in predicting the 
evolution of a disease in a dynamic contact network. In 
future work, this model can be adapted and extended to 
account for individuals cutting and creating links with 
knowledge of the state of others in the population, i.e., 
link-type dependent network dynamics. This two-way 
feedback will lead to more sophisticated network proper- 
ties such as degree correlations, high clustering or even 
network fragmentation. In such cases ODE models need 
to be used with care, making sure that the agreement 
with simulations remains valid. Besides modelling epi- 
demics, this framework could also be used to study the 
spread of information, beliefs and new ideas within popu- 
lations, and as such could have implications across a wide 
range of disciplines beyond the mathematical biology and 
physics communities. 
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